clear all;

num_samples = 20;
%num_samples = 2;
cvx_quiet(0)
cvx_solver sdpt3
tol = 1e-2;

gammas = repmat(0:0.005:0.1, 1, num_samples);
%gammas = repmat(0:0.05:0.1, 1, num_samples);
Bregs = cell(2, length(gammas));
Gregs = cell(2, length(gammas));
aic_regs = zeros(2, length(gammas));
bic_regs = zeros(2, length(gammas));
sens = zeros(2, length(gammas));
spec = zeros(2, length(gammas));
norm_errs = zeros(2, length(gammas));
for i=1:length(gammas)
	if mod(i,length(unique(gammas))) == 1
		i_old = i;
		generate;
		i = i_old;
		Bbase = full(B);
		Ybase = compute_D(B'*B,n,p);
		Gbase = ones(n);
		for j=0:p
			Gbase = Gbase & abs(Ybase(:,(1:n) + n*j) < tol);
		end
		Gbase = ~Gbase;
	end

	gamma = gammas(i);

	Bmlm = regular_mlm(C, n, p, gamma); 
	Ymlm = compute_D(Bmlm'*Bmlm, n, p); 
	Gmlm = ones(n);
	for j=0:p
		Gmlm = Gmlm & abs(Ymlm(:,(1:n) + n*j) < tol);
	end
	Bregs{1,i} = Bmlm;
	Gregs{1,i} = ~Gmlm;
	aic_regs(1,i) = aic(Bmlm'*Bmlm, C, n, p, N, tol);
	bic_regs(1,i) = bic(Bmlm'*Bmlm, C, n, p, N, tol);
	sens(1,i) = compute_sensitivity(Ybase, Ymlm, tol);
	spec(1,i) = compute_specificity(Ybase, Ymlm, tol);
	norm_errs(1,i) = norm(Bmlm-Bbase, 'fro');
	fprintf('computed coefficients by maximum likelihood; gamma = %1.3f\n', gamma);

%	Bhyb = regular_hyb(C, n, p, gamma, tol, 1e3); 
%	Yhyb = compute_D(Bhyb'*Bhyb, n, p);
%	aic_regs(2,i) = aic(Bhyb'*Bhyb, C, n, p, N, tol);
%	bic_regs(2,i) = bic(Bhyb'*Bhyb, C, n, p, N, tol);
%	sens(2,i) = compute_sensitivity(Ybase, Yhyb, tol);
%	spec(2,i) = compute_specificity(Ybase, Yhyb, tol);
%	fprintf('Computed coefficients by regression on the AR coefficients and ML on Sigma\n')
%
	Brgr = regular_rgr(C, n, p, gamma, tol, 1e3); 
	Yrgr = compute_D(Brgr'*Brgr, n, p);
	Grgr = ones(n);
	for j=0:p
		Grgr = Grgr & abs(Yrgr(:,(1:n) + n*j) < tol);
	end
	Bregs{2,i} = Brgr;
	Gregs{2,i} = ~Grgr;
	aic_regs(2,i) = aic(Brgr'*Brgr, C, n, p, N, tol);
	bic_regs(2,i) = bic(Brgr'*Brgr, C, n, p, N, tol);
	sens(2,i) = compute_sensitivity(Ybase, Yrgr, tol);
	spec(2,i) = compute_specificity(Ybase, Yrgr, tol);
	norm_errs(2,i) = norm(Brgr-Bbase, 'fro');
	fprintf('Computed coefficients by iterative regression; gamma = %1.3f\n', gamma)

%	%tic;
%	Bdul = dual_reg_rgr(C, n, p, gamma); 
%	Ydul = compute_D(Bdul'*Bdul, n, p);
%	aic_regs(4,i) = aic(Bdul'*Bdul, C, n, p, N, tol);
%	bic_regs(4,i) = bic(Bdul'*Bdul, C, n, p, N, tol);
%	sens(4,i) = compute_sensitivity(Ybase, Ydul, tol);
%	spec(4,i) = compute_specificity(Ybase, Ydul, tol);
%	%tdul=toc; 
%	fprintf('Computed coefficients by dual of regression\n')
end
%for i=1:2                                                           
%	for j=1:21                                                          
%		dlmwrite(sprintf('../Data/graph_n%i_g%1.3f.csv', n, gammas(j)), Gregs{i,j})
%	end
%end

p_max = 6; % TODO this is arbitrary
C_max = windowed_est(x, p_max);
lambdas = 0.05:0.01:0.2;
Bblocks = cell(size(lambdas));
p_opts = zeros(size(lambdas));
aic_blocks = zeros(size(lambdas));
bic_blocks = zeros(size(lambdas));
%[bic_min, i] = min(mean(bic_regs));
%gamma = gammas(i);
gamma = 0.015;
for i=1:length(lambdas)
	lambda_0 = lambdas(i);
	Bblock = regular_block_mlm(C_max, n, p_max, gamma, lambda_0);
	Yblock = compute_D(Bblock'*Bblock, n, p_max);
	[nz_rows,nz_cols] = find(abs(Yblock) > tol);
	p_opt = ceil(max(nz_cols) / n) - 1;
	Bblock = Bblock(1:n,1:(n*(p_opt+1)));
	Yblock = compute_D(Bblock'*Bblock, n, p_opt);
	Xblock = Bblock'*Bblock;
	C_opt = C_max(1:(n*(p_opt+1)),1:(n*(p_opt+1)));
%	Yblocked = blockify(Yblock);
%	Gblock = abs(Yblock) > tol;
	Bblocks{i} = Bblock;
	p_opts(i) = p_opt;
	aic_blocks(i) = aic(Xblock, C_opt, n, p_opt, N, tol);
	bic_blocks(i) = bic(Xblock, C_opt, n, p_opt, N, tol);
	fprintf('Computed coefficients with block peanlty; gamma=%1.3f, lambda=%1.3f\n', gamma, lambda_0)
end
save(sprintf('../Data/n%i_p%i_N%i.mat', n, p, N));
[spec_smooth_1, sens_smooth_1] = smooth_values(spec(1,:),sens(1,:));
[spec_smooth_2, sens_smooth_2] = smooth_values(spec(2,:),sens(2,:));
spec_smooth = [spec_smooth_1; spec_smooth_2];
sens_smooth = [sens_smooth_1; sens_smooth_2];
%plot(1-spec_smooth',sens_smooth')
%plot(gammas, norm_errs'u)

%figure;
%[ax, h1, h2] = plotyy(lambdas, aic_blocks, lambdas, p_opts);
%default_lims = axis(ax(2));
%axis(ax(2), [default_lims(1:2) min(p_opts)-0.5 max(p_opts)+0.5]);
%hold on
%plot(lambdas, bic_blocks, 'r');
%axis auto
%legend('AIC','BIC','K_{est}', 'Location', 'Best')
%legend boxoff  
